Metrics for assessing community assembly processes

source('../env.R')

community_data = read_csv(filename(COMMUNITY_OUTPUT_DIR, 'community_assembly_metrics.csv'))

Rows: 792 Columns: 93── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr  (1): is_urban_threshold
dbl (92): city_id, mntd_normalised, mntd_actual, mntd_min, mntd_max, mntd_mean, mntd_sd, fd_normalised, fd_actual, fd_min, fd_max, fd_mean, fd_sd, FRic_normal...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

community_data$is_urban_threshold = factor(community_data$is_urban_threshold, levels = c('low', 'medium', 'high'), labels = c('Low', 'Medium', 'High'))
head(community_data)

colnames(community_data)

 [1] "city_id"                           "mntd_normalised"                   "mntd_actual"                       "mntd_min"                         
 [5] "mntd_max"                          "mntd_mean"                         "mntd_sd"                           "fd_normalised"                    
 [9] "fd_actual"                         "fd_min"                            "fd_max"                            "fd_mean"                          
[13] "fd_sd"                             "FRic_normalised"                   "FRic_actual"                       "FRic_min"                         
[17] "FRic_max"                          "FRic_mean"                         "FRic_sd"                           "mass_fd_normalised"               
[21] "mass_fd_actual"                    "mass_fd_min"                       "mass_fd_max"                       "mass_fd_mean"                     
[25] "mass_fd_sd"                        "mass_var_normalised"               "mass_var_actual"                   "mass_var_min"                     
[29] "mass_var_max"                      "mass_var_mean"                     "mass_var_sd"                       "mass_sdndr_normalised"            
[33] "mass_sdndr_actual"                 "mass_sdndr_min"                    "mass_sdndr_max"                    "mass_sdndr_mean"                  
[37] "mass_sdndr_sd"                     "locomotory_trait_fd_normalised"    "locomotory_trait_fd_actual"        "locomotory_trait_fd_min"          
[41] "locomotory_trait_fd_max"           "locomotory_trait_fd_mean"          "locomotory_trait_fd_sd"            "locomotory_trait_var_normalised"  
[45] "locomotory_trait_var_actual"       "locomotory_trait_var_min"          "locomotory_trait_var_max"          "locomotory_trait_var_mean"        
[49] "locomotory_trait_var_sd"           "locomotory_trait_sdndr_normalised" "locomotory_trait_sdndr_actual"     "locomotory_trait_sdndr_min"       
[53] "locomotory_trait_sdndr_max"        "locomotory_trait_sdndr_mean"       "locomotory_trait_sdndr_sd"         "trophic_trait_fd_normalised"      
[57] "trophic_trait_fd_actual"           "trophic_trait_fd_min"              "trophic_trait_fd_max"              "trophic_trait_fd_mean"            
[61] "trophic_trait_fd_sd"               "trophic_trait_var_normalised"      "trophic_trait_var_actual"          "trophic_trait_var_min"            
[65] "trophic_trait_var_max"             "trophic_trait_var_mean"            "trophic_trait_var_sd"              "trophic_trait_sdndr_normalised"   
[69] "trophic_trait_sdndr_actual"        "trophic_trait_sdndr_min"           "trophic_trait_sdndr_max"           "trophic_trait_sdndr_mean"         
[73] "trophic_trait_sdndr_sd"            "gape_width_fd_normalised"          "gape_width_fd_actual"              "gape_width_fd_min"                
[77] "gape_width_fd_max"                 "gape_width_fd_mean"                "gape_width_fd_sd"                  "gape_width_var_normalised"        
[81] "gape_width_var_actual"             "gape_width_var_min"                "gape_width_var_max"                "gape_width_var_mean"              
[85] "gape_width_var_sd"                 "gape_width_sdndr_normalised"       "gape_width_sdndr_actual"           "gape_width_sdndr_min"             
[89] "gape_width_sdndr_max"              "gape_width_sdndr_mean"             "gape_width_sdndr_sd"               "urban_pool_size"                  
[93] "is_urban_threshold"               

Join on realms

city_to_realm = read_csv(filename(CITY_DATA_OUTPUT_DIR, 'realms.csv'))

Rows: 342 Columns: 2── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (1): core_realm
dbl (1): city_id
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

community_data_with_realm = left_join(community_data, city_to_realm)

Joining with `by = join_by(city_id)`

Cities as points

city_points = st_centroid(read_sf(filename(CITY_DATA_OUTPUT_DIR, 'city_selection.shp')))

Warning: st_centroid assumes attributes are constant over geometries of x

city_points_low = city_points %>% left_join(community_data[community_data$is_urban_threshold == 'Low',])

Joining with `by = join_by(city_id)`

city_points_med = city_points %>% left_join(community_data[community_data$is_urban_threshold == 'Medium',])

Joining with `by = join_by(city_id)`

city_points_high = city_points %>% left_join(community_data[community_data$is_urban_threshold == 'High',])

Joining with `by = join_by(city_id)`

sf::sf_use_s2(FALSE)

Spherical geometry (s2) switched off

COUNTRY_BOUNDARIES = '/Users/james/Dropbox/PhD/WorldBank_countries_Admin0_10m/WB_countries_Admin0_10m.shp'
world_map = st_simplify(st_read(COUNTRY_BOUNDARIES), dTolerance = 0.02)

Reading layer `WB_countries_Admin0_10m' from data source `/Users/james/Dropbox/PhD/WorldBank_countries_Admin0_10m/WB_countries_Admin0_10m.shp' using driver `ESRI Shapefile'
Simple feature collection with 251 features and 52 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -180 ymin: -59.47275 xmax: 180 ymax: 83.6341
Geodetic CRS:  WGS 84

Warning: st_simplify does not correctly simplify longitude/latitude data, dTolerance needs to be in decimal degrees

normalised_colours_scale = scale_colour_gradient2(
    low = "darkgreen",
    mid = "yellow",
    high = "red",
    midpoint = 0.5,
    space = "Lab",
    na.value = "grey50",
    guide = "colourbar",
    aesthetics = "colour",
  )

Load community data, and create long format version

communities = read_csv(filename(COMMUNITY_OUTPUT_DIR, 'communities_for_analysis.csv'))

Rows: 2462 Columns: 9── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (5): city_name, jetz_species_name, seasonal, presence, origin
dbl (1): city_id
lgl (3): present_urban_high, present_urban_med, present_urban_low
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

communities_summary = communities %>% group_by(city_id) %>% summarise(
  regional_pool_size = n(), 
  urban_size_high = sum(present_urban_high), 
  urban_size_med = sum(present_urban_med), 
  urban_size_low = sum(present_urban_low)
)

communities_summary_long = bind_rows(
  communities_summary %>% rename(urban_pool_size = 'urban_size_high') %>% dplyr::select(city_id, regional_pool_size, urban_pool_size) %>% mutate(is_urban_threshold = 'High'),
  communities_summary %>% rename(urban_pool_size = 'urban_size_med') %>% dplyr::select(city_id, regional_pool_size, urban_pool_size) %>% mutate(is_urban_threshold = 'Medium'),
  communities_summary %>% rename(urban_pool_size = 'urban_size_low') %>% dplyr::select(city_id, regional_pool_size, urban_pool_size) %>% mutate(is_urban_threshold = 'Low')
)
communities_summary_long

Load trait data

traits = read_csv(filename(TAXONOMY_OUTPUT_DIR, 'traits_jetz.csv'))

Rows: 304 Columns: 5── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (1): jetz_species_name
dbl (4): gape_width, trophic_trait, locomotory_trait, mass
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

head(traits)

Load realm geo

st_crs(resolve)

Coordinate Reference System:
  User input: WGS 84 
  wkt:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["latitude",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["longitude",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]

Examine individual metrics

geom_normalised_histogram = function(name, gg, legend.position = "right") {
  gg + 
    geom_histogram(aes(fill = core_realm), binwidth = 0.1, position = "dodge") +
    geom_vline(aes(xintercept = 0.5), color = "#000000", size = 0.4) +
    geom_vline(aes(xintercept = 0), color = "#000000", size = 0.2, linetype = "dashed") +
    geom_vline(aes(xintercept = 1), color = "#000000", size = 0.2, linetype = "dashed") + 
    ylab("Number of cities") + xlab("Normalised Response") + ylim(c(0, 70)) +
    labs(title = name, fill = 'Realm') +
    facet_wrap(~ is_urban_threshold, ncol=1) + 
    theme_bw() +
    theme(legend.position=legend.position)
}

geom_map = function(map_sf, title) {
  norm_mntd_analysis_geo = ggplot() + 
    geom_sf(data = world_map, aes(geometry = geometry)) +
    map_sf +
    normalised_colours_scale +
    labs(title = title, colour = 'Normalised\nResponse')
}

MNTD

norm_mntd_analysis_plot = geom_normalised_histogram(
  'MNTD', 
  ggplot(community_data_with_realm, aes(mntd_normalised))
)

Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
Please use `linewidth` instead.

norm_mntd_analysis_plot

norm_mntd_analysis_geo = geom_map(geom_sf(data = city_points_low, aes(color = mntd_normalised, geometry = geometry)), 'Low threshold')
norm_mntd_analysis_geo

ggarrange(norm_mntd_analysis_plot, 
          ggarrange(
            geom_normalised_histogram('MNTD', ggplot(community_data_with_realm, aes(mntd_normalised)), legend.position = "bottom"), 
            geom_map(geom_sf(data = city_points_med, aes(color = mntd_normalised, geometry = geometry)), 'Medium threshold'),
            geom_map(geom_sf(data = city_points_high, aes(color = mntd_normalised, geometry = geometry)), 'High threshold'),
            ncol=1, common.legend = T),
          ncol = 2)

ggsave(filename(FIGURES_OUTPUT_DIR, 'normalised_mntd.jpg'), width = 2500, height=2500, units = 'px')

FD

norm_fd_analysis_plot = geom_normalised_histogram(
  'FD', 
  ggplot(community_data_with_realm, aes(fd_normalised))
)
norm_fd_analysis_plot

norm_fd_analysis_geo = geom_map(geom_sf(data = city_points_low, aes(color = fd_normalised, geometry = geometry)), 'Low threshold')
norm_fd_analysis_geo

ggarrange(geom_normalised_histogram('FD', ggplot(community_data_with_realm, aes(fd_normalised)), legend.position = "bottom"), 
          ggarrange(
            norm_fd_analysis_geo, 
            geom_map(geom_sf(data = city_points_med, aes(color = fd_normalised, geometry = geometry)), 'Medium threshold'),
            geom_map(geom_sf(data = city_points_high, aes(color = fd_normalised, geometry = geometry)), 'High threshold'),
            ncol=1, common.legend = T),
          ncol = 2)

ggsave(filename(FIGURES_OUTPUT_DIR, 'normalised_fd.jpg'), width = 2500, height=2500, units = 'px')

FD - Locomotory

norm_fd_loco_analysis_plot = geom_normalised_histogram(
  'FD Locomotory', 
  ggplot(community_data_with_realm, aes(locomotory_trait_fd_normalised))
)
norm_fd_loco_analysis_plot

VAR - Locomotory

norm_var_loco_analysis_plot = geom_normalised_histogram(
  'Var Locomotory', 
  ggplot(community_data_with_realm, aes(locomotory_trait_var_normalised))
)
norm_var_loco_analysis_plot

FD - Trophic

norm_fd_trophic_analysis_plot = geom_normalised_histogram(
  'FD Trophic', 
  ggplot(community_data_with_realm, aes(trophic_trait_fd_normalised))
)
norm_fd_trophic_analysis_plot

VAR - Trophic

norm_var_trophic_analysis_plot = geom_normalised_histogram(
  'Var Trophic', 
  ggplot(community_data_with_realm, aes(trophic_trait_var_normalised))
)
norm_var_trophic_analysis_plot

FD - Gape Width

norm_fd_gape_analysis_plot = geom_normalised_histogram(
  'FD Gape Width', 
  ggplot(community_data_with_realm, aes(gape_width_fd_normalised))
)
norm_fd_gape_analysis_plot

VAR - Gape Width

norm_var_gape_analysis_plot = geom_normalised_histogram(
  'Var Gape Width', 
  ggplot(community_data_with_realm, aes(gape_width_var_normalised))
)
norm_var_gape_analysis_plot

FD - Mass

norm_fd_mass_analysis_plot = geom_normalised_histogram(
  'FD Mass', 
  ggplot(community_data_with_realm, aes(mass_fd_normalised))
)
norm_fd_mass_analysis_plot

VAR - Mass

norm_var_mass_analysis_plot = geom_normalised_histogram(
  'Var Mass', 
  ggplot(community_data_with_realm, aes(mass_var_normalised))
)
norm_var_mass_analysis_plot

Compare metrics against each other

ggplot(community_data_with_realm, aes(x = fd_normalised, y = mntd_normalised, colour = core_realm)) + 
  geom_point() +
  ylab("MNTD") + 
  xlab("FD") +
  theme_bw() +
  facet_wrap(~ is_urban_threshold, ncol = 2) + labs(color = "Realm")
ggsave(filename(FIGURES_OUTPUT_DIR, 'mntd_by_fd.jpg'))

Saving 7.29 x 4.51 in image

mntd_fd_analysis = community_data_with_realm %>% 
  dplyr::select(city_id, is_urban_threshold, mntd_normalised, fd_normalised) %>%
  left_join(communities_summary_long) %>%
  mutate(urban_pool_perc = urban_pool_size * 100 / regional_pool_size)

Joining with `by = join_by(city_id, is_urban_threshold)`

mntd_fd_analysis

ggpairs(mntd_fd_analysis %>% dplyr::filter(is_urban_threshold == 'Low') %>% dplyr::select(mntd_normalised, fd_normalised, regional_pool_size, urban_pool_size, urban_pool_perc))

ggpairs(mntd_fd_analysis %>% dplyr::filter(is_urban_threshold == 'Medium') %>% dplyr::select(mntd_normalised, fd_normalised, regional_pool_size, urban_pool_size, urban_pool_perc))
ggsave(filename(FIGURES_OUTPUT_DIR, 'mntd_by_fd_pairs_medium.jpg'))

Saving 7.29 x 4.51 in image

ggpairs(mntd_fd_analysis %>% dplyr::filter(is_urban_threshold == 'High') %>% dplyr::select(mntd_normalised, fd_normalised, regional_pool_size, urban_pool_size, urban_pool_perc))

Plot Urban/Regional Traits

communities_with_traits = communities %>% left_join(traits) %>% left_join(city_to_realm) %>% filter(core_realm != 'Oceania')

Joining with `by = join_by(jetz_species_name)`Joining with `by = join_by(city_id)`

head(communities_with_traits)

realms = unique(communities_with_traits$core_realm)
realms

[1] "Nearctic"    "Neotropic"   "Palearctic"  "Afrotropic"  "Indomalayan" "Australasia"

communities_with_traits_long = bind_rows(
  communities_with_traits %>% dplyr::select(gape_width, trophic_trait, locomotory_trait, mass, core_realm) %>% mutate(record_type = 'Regional'),
  communities_with_traits %>% filter(present_urban_low) %>% dplyr::select(gape_width, trophic_trait, locomotory_trait, mass, core_realm) %>% mutate(record_type = 'Urban Low Threshold'),
  communities_with_traits %>% filter(present_urban_med) %>% dplyr::select(gape_width, trophic_trait, locomotory_trait, mass, core_realm) %>% mutate(record_type = 'Urban Medium Threshold'),
  communities_with_traits %>% filter(present_urban_high) %>% dplyr::select(gape_width, trophic_trait, locomotory_trait, mass, core_realm) %>% mutate(record_type = 'Urban High Threshold')
)
communities_with_traits_long$record_type = factor(communities_with_traits_long$record_type, levels = c('Regional', 'Urban Low Threshold', 'Urban Medium Threshold', 'Urban High Threshold'))

head(communities_with_traits_long)

convex_hull_loco_trophic_per_realm = function(filtered_df, realms) {
  result = data.frame()
  
  for (realm in realms) {
    result = rbind(result, 
      filtered_df %>% 
        filter(core_realm == realm) %>% 
        slice(chull(trophic_trait, locomotory_trait)) %>% 
        dplyr::select(trophic_trait, locomotory_trait) %>%
        mutate(core_realm = realm)
    )
  }
  
  result
}

regional_hull = convex_hull_loco_trophic_per_realm(communities_with_traits, realms)
urban_hull_high = convex_hull_loco_trophic_per_realm(communities_with_traits %>% filter(present_urban_high), realms)
urban_hull_med = convex_hull_loco_trophic_per_realm(communities_with_traits %>% filter(present_urban_med), realms)
urban_hull_low = convex_hull_loco_trophic_per_realm(communities_with_traits %>% filter(present_urban_low), realms)

ggplot(data = communities_with_traits_long, aes(x = trophic_trait, y = locomotory_trait)) +
  geom_polygon(data = regional_hull, alpha = 0.1, fill = "green", color="green") +
  geom_polygon(data = urban_hull_low, alpha = 0.15, fill = "yellow", color="yellow") +
  geom_polygon(data = urban_hull_med, alpha = 0.2, fill = "orange", color = "orange") +
  geom_polygon(data = urban_hull_high, alpha = 0.25, fill = "red", color = "red") +
  geom_point(aes(colour = record_type)) +
  theme_bw() + scale_color_manual(values=c("green", "yellow", "orange", "red")) +
  xlab('Trophic Trait') + ylab('Locomotory Trait') + 
  theme(legend.position="bottom") + 
  facet_wrap(~ core_realm)
ggsave(filename(FIGURES_OUTPUT_DIR, 'traits_by_realm_loco_trophic.jpg'))

Saving 7.29 x 4.51 in image

convex_hull_gape_mass_per_realm = function(filtered_df, realms) {
  result = data.frame()
  
  for (realm in realms) {
    result = rbind(result, 
      filtered_df %>% 
        filter(core_realm == realm) %>% 
        slice(chull(gape_width, mass)) %>% 
        dplyr::select(gape_width, mass) %>%
        mutate(core_realm = realm)
    )
  }
  
  result
}

regional_hull_gm = convex_hull_gape_mass_per_realm(communities_with_traits, realms)
urban_hull_high_gm = convex_hull_gape_mass_per_realm(communities_with_traits %>% filter(present_urban_high), realms)
urban_hull_med_gm = convex_hull_gape_mass_per_realm(communities_with_traits %>% filter(present_urban_med), realms)
urban_hull_low_gm = convex_hull_gape_mass_per_realm(communities_with_traits %>% filter(present_urban_low), realms)

ggplot(data = communities_with_traits_long, aes(x = gape_width, y = mass)) +
  geom_polygon(data = regional_hull_gm, alpha = 0.1, fill = "green", color="green") +
  geom_polygon(data = urban_hull_low_gm, alpha = 0.15, fill = "yellow", color="yellow") +
  geom_polygon(data = urban_hull_med_gm, alpha = 0.2, fill = "orange", color = "orange") +
  geom_polygon(data = urban_hull_high_gm, alpha = 0.25, fill = "red", color = "red") +
  geom_point(aes(colour = record_type)) +
  theme_bw() + scale_color_manual(values=c("green", "yellow", "orange", "red")) +
  xlab('Gape Width') + ylab('Mass') + 
  theme(legend.position="bottom") + 
  facet_wrap(~ core_realm)
ggsave(filename(FIGURES_OUTPUT_DIR, 'traits_by_realm_gape_mass.jpg'))

Saving 7.29 x 4.51 in image

Plot Urban/Regional Phylogeny

phylo_tree = read.tree(filename(TAXONOMY_OUTPUT_DIR, 'phylogeny.tre'))
ggtree(phylo_tree, layout='circular')

plot_phylo = function(realm, community_df = communities_with_traits) {
  species_df = community_df %>% filter(core_realm == realm) %>% group_by(jetz_species_name) %>% summarise(
    regional_pools = n(),
    urban_pools_high = sum(present_urban_high),
    urban_pools_medium = sum(present_urban_med),
    urban_pools_low = sum(present_urban_low)
  ) %>% mutate(
    urban_tolerence_high = urban_pools_high / regional_pools,
    urban_tolerence_med = urban_pools_medium / regional_pools,
    urban_tolerence_low = urban_pools_low / regional_pools
  )
  
  tree_cropped <- ladderize(drop.tip(phylo_tree, setdiff(phylo_tree$tip.label, species_df$jetz_species_name)))
  
  p = ggtree(tree_cropped) + geom_tiplab(align=T) 
  p1 = facet_plot(p, panel='High', data=species_df, geom=geom_segment, aes(x=0, xend=urban_tolerence_high, y=y, yend=y), size=3, color='red') 
  p2 = facet_plot(p1, panel='Medium', data=species_df, geom=geom_segment, aes(x=0, xend=urban_tolerence_med, y=y, yend=y), size=3, color='orange') 
  p3 = facet_plot(p2, panel='Low', data=species_df, geom=geom_segment, aes(x=0, xend=urban_tolerence_low, y=y, yend=y), size=3, color='yellow') 
  
  facet_widths(p3 + xlim_tree(60) + theme_tree2(), c(Tree = 5)) + labs(title = realm, subtitle = 'Urban tolerance')
}

realms

[1] "Nearctic"    "Neotropic"   "Palearctic"  "Afrotropic"  "Indomalayan" "Australasia"

plot_phylo('Nearctic')

ggsave(filename(FIGURES_OUTPUT_DIR, 'phylogeny_nearctic.jpg'))

Saving 7.29 x 4.51 in image

plot_phylo('Neotropic')

ggsave(filename(FIGURES_OUTPUT_DIR, 'phylogeny_neotropic.jpg'), width=2187, units="px")

Saving 2187 x 3000 px image

plot_phylo('Palearctic')

ggsave(filename(FIGURES_OUTPUT_DIR, 'phylogeny_palearctic.jpg'))

Saving 7.29 x 4.51 in image

plot_phylo('Afrotropic')

ggsave(filename(FIGURES_OUTPUT_DIR, 'phylogeny_afrotropic.jpg'))

Saving 7.29 x 4.51 in image

plot_phylo('Australasia')

ggsave(filename(FIGURES_OUTPUT_DIR, 'phylogeny_australasia.jpg'))

Saving 7.29 x 4.51 in image

Species in communities

It seems reasonable to expect that cities with simialr regional pools will have similar species entering the city, and thus a similar response to urbanisation.

to_species_matrix = function(filtered_communities) {
  filtered_communities %>% 
    dplyr::select(city_id, jetz_species_name) %>% 
    distinct() %>%
    mutate(present = TRUE) %>% 
    pivot_wider(
      names_from = jetz_species_name, 
      values_from = "present", 
      values_fill = list(present = F)
    ) %>% 
    tibble::column_to_rownames(var='city_id')
}

community_nmds = function(filtered_communities) {
  species_matrix = to_species_matrix(filtered_communities)
  nmds <- metaMDS(species_matrix, k=2, trymax = 30)
  nmds_result = data.frame(scores(nmds))
  nmds_result$city_id = as.double(rownames(scores(nmds)))
  rownames(nmds_result) = NULL
  
  nmds_result %>%
    left_join(city_points[,c('city_id', 'city_nm')])
}

plot_nmds = function(plot_data, clusters, realm, plot_is_urban_threshold) {
  pd = plot_data %>%
    left_join(community_data) %>%
    filter(is_urban_threshold == plot_is_urban_threshold)
  
  ggplot(pd, aes(x = NMDS1, y = NMDS2)) + 
    geom_polygon(data = clusters, aes(fill = cluster), alpha = 0.2) +
    geom_jitter(aes(colour = mntd_normalised, size = fd_normalised)) +
    normalised_colours_scale +
    geom_text_repel(aes(label = city_nm), max.overlaps = 25) +
    labs(
      title = realm, 
      subtitle = paste('Regional pool, with urban community values using', plot_is_urban_threshold, 'threshold'), 
      color = 'MNTD Normalised', 
      size = 'FD Normalised')
}

https://www.datacamp.com/tutorial/k-means-clustering-r

scree_plot = function(community_nmds_data) {
  # Decide how many clusters to look at
  n_clusters <- 10
  
  # Initialize total within sum of squares error: wss
  wss <- numeric(n_clusters)
  
  set.seed(123)
  
  # Look over 1 to n possible clusters
  for (i in 1:n_clusters) {
    # Fit the model: km.out
    km.out <- kmeans(community_nmds_data[,c('NMDS1','NMDS2')], centers = i, nstart = 20)
    # Save the within cluster sum of squares
    wss[i] <- km.out$tot.withinss
  }
  
  # Produce a scree plot
  wss_df <- tibble(clusters = 1:n_clusters, wss = wss)
   
  scree_plot <- ggplot(wss_df, aes(x = clusters, y = wss, group = 1)) +
      geom_point(size = 4)+
      geom_line() +
      scale_x_continuous(breaks = c(2, 4, 6, 8, 10)) +
      xlab('Number of clusters')
  scree_plot
}

convex_hull_by_cluster = function(community_df) {
  result = data.frame()
  
  clusters = unique(community_df$cluster)
  for (cluster_i in clusters) {
    result = rbind(result, 
      community_df %>% 
        filter(cluster == cluster_i) %>% 
        slice(chull(NMDS1, NMDS2)) %>% 
        dplyr::select(NMDS1, NMDS2, cluster) 
    )
  }
  
  result
}

Nearctic

nearctic_communities = community_nmds(communities_with_traits %>% filter(core_realm == 'Nearctic'))    

Run 0 stress 0.06772313 
Run 1 stress 0.09925753 
Run 2 stress 0.0814111 
Run 3 stress 0.0714024 
Run 4 stress 0.08477133 
Run 5 stress 0.06792384 
... Procrustes: rmse 0.002677399  max resid 0.01592546 
Run 6 stress 0.07140237 
Run 7 stress 0.09410676 
Run 8 stress 0.06772314 
... Procrustes: rmse 0.000007756193  max resid 0.00001712457 
... Similar to previous best
Run 9 stress 0.0847715 
Run 10 stress 0.06772321 
... Procrustes: rmse 0.00005153412  max resid 0.00008951089 
... Similar to previous best
Run 11 stress 0.06792388 
... Procrustes: rmse 0.00267859  max resid 0.01598036 
Run 12 stress 0.0677232 
... Procrustes: rmse 0.00004962282  max resid 0.00008595323 
... Similar to previous best
Run 13 stress 0.08193052 
Run 14 stress 0.09279085 
Run 15 stress 0.08494123 
Run 16 stress 0.09422939 
Run 17 stress 0.07122361 
Run 18 stress 0.08521544 
Run 19 stress 0.07427005 
Run 20 stress 0.06792383 
... Procrustes: rmse 0.002676712  max resid 0.01589483 
*** Solution reached

Joining with `by = join_by(city_id)`

nearctic_communities

scree_plot(nearctic_communities)

Warning: Quick-TRANSfer stage steps exceeded maximum (= 3300)Warning: Quick-TRANSfer stage steps exceeded maximum (= 3300)Warning: Quick-TRANSfer stage steps exceeded maximum (= 3300)Warning: Quick-TRANSfer stage steps exceeded maximum (= 3300)Warning: did not converge in 10 iterationsWarning: Quick-TRANSfer stage steps exceeded maximum (= 3300)Warning: Quick-TRANSfer stage steps exceeded maximum (= 3300)Warning: did not converge in 10 iterationsWarning: Quick-TRANSfer stage steps exceeded maximum (= 3300)Warning: Quick-TRANSfer stage steps exceeded maximum (= 3300)Warning: did not converge in 10 iterationsWarning: Quick-TRANSfer stage steps exceeded maximum (= 3300)Warning: did not converge in 10 iterationsWarning: Quick-TRANSfer stage steps exceeded maximum (= 3300)Warning: did not converge in 10 iterationsWarning: Quick-TRANSfer stage steps exceeded maximum (= 3300)Warning: Quick-TRANSfer stage steps exceeded maximum (= 3300)Warning: did not converge in 10 iterationsWarning: Quick-TRANSfer stage steps exceeded maximum (= 3300)Warning: Quick-TRANSfer stage steps exceeded maximum (= 3300)Warning: Quick-TRANSfer stage steps exceeded maximum (= 3300)Warning: did not converge in 10 iterationsWarning: Quick-TRANSfer stage steps exceeded maximum (= 3300)

nearctic_clusters <- kmeans(nearctic_communities[,c('NMDS1', 'NMDS2')], centers = 3, nstart = 20)
nearctic_communities$cluster = as.factor(nearctic_clusters$cluster)

ggplot(nearctic_communities, aes(x = NMDS1, y = NMDS2, colour = cluster)) + geom_point()

nearctic_cluster_hulls = convex_hull_by_cluster(nearctic_communities)

plot_nmds(nearctic_communities, nearctic_cluster_hulls, 'Nearctic', 'High')

Joining with `by = join_by(city_id)`

plot_nmds(nearctic_communities, nearctic_cluster_hulls, 'Nearctic', 'Medium')

Joining with `by = join_by(city_id)`

plot_nmds(nearctic_communities, nearctic_cluster_hulls, 'Nearctic', 'Low')

Joining with `by = join_by(city_id)`

species_in_neartic = nearctic_communities %>% dplyr::select(city_id, cluster) %>% left_join(communities_with_traits)

Joining with `by = join_by(city_id)`

plot_phylo('Nearctic', species_in_neartic %>% filter(cluster == 1))

ggplot() +
  geom_sf(data = resolve_nearctic, aes(geometry = geometry), alpha = 0.2) +
  geom_sf(data = nearctic_communities %>% filter(cluster == 1) %>% left_join(city_points_med), aes(color = mntd_normalised, geometry = geometry)) +
  normalised_colours_scale + theme_bw()

Joining with `by = join_by(city_id, city_nm, geometry)`

plot_phylo('Nearctic', species_in_neartic %>% filter(cluster == 2))

ggplot() +
  geom_sf(data = resolve_nearctic, aes(geometry = geometry), alpha = 0.2) +
  geom_sf(data = nearctic_communities %>% filter(cluster == 2) %>% left_join(city_points_med), aes(color = mntd_normalised, geometry = geometry)) +
  normalised_colours_scale + theme_bw()

Joining with `by = join_by(city_id, city_nm, geometry)`

ggplot() +
  geom_sf(data = resolve_nearctic, aes(geometry = geometry), alpha = 0.2) +
  geom_sf(data = nearctic_communities %>% filter(cluster == 3) %>% left_join(city_points_med), aes(color = mntd_normalised, geometry = geometry)) +
  normalised_colours_scale + theme_bw()

Joining with `by = join_by(city_id, city_nm, geometry)`

Neotropic

neotropic_communities = community_nmds(communities_with_traits %>% filter(core_realm == 'Neotropic'))    

Run 0 stress 0.134619 
Run 1 stress 0.1344509 
... New best solution
... Procrustes: rmse 0.006928659  max resid 0.0457892 
Run 2 stress 0.1369387 
Run 3 stress 0.134619 
... Procrustes: rmse 0.006940788  max resid 0.04575861 
Run 4 stress 0.1344509 
... New best solution
... Procrustes: rmse 0.00004403156  max resid 0.0001387448 
... Similar to previous best
Run 5 stress 0.1346406 
... Procrustes: rmse 0.007192008  max resid 0.0457332 
Run 6 stress 0.1346226 
... Procrustes: rmse 0.007280471  max resid 0.04573126 
Run 7 stress 0.136377 
Run 8 stress 0.1348046 
... Procrustes: rmse 0.006663518  max resid 0.04610571 
Run 9 stress 0.1406945 
Run 10 stress 0.1346226 
... Procrustes: rmse 0.007291867  max resid 0.04571052 
Run 11 stress 0.134619 
... Procrustes: rmse 0.006932652  max resid 0.04571971 
Run 12 stress 0.134433 
... New best solution
... Procrustes: rmse 0.001186888  max resid 0.006991429 
... Similar to previous best
Run 13 stress 0.1348237 
... Procrustes: rmse 0.006651112  max resid 0.04607972 
Run 14 stress 0.1346405 
... Procrustes: rmse 0.007243159  max resid 0.04569399 
Run 15 stress 0.1406945 
Run 16 stress 0.134619 
... Procrustes: rmse 0.006830928  max resid 0.04574414 
Run 17 stress 0.1344509 
... Procrustes: rmse 0.001186553  max resid 0.00692233 
... Similar to previous best
Run 18 stress 0.141222 
Run 19 stress 0.1346207 
... Procrustes: rmse 0.006999734  max resid 0.04555 
Run 20 stress 0.134433 
... New best solution
... Procrustes: rmse 0.00002505145  max resid 0.00008159977 
... Similar to previous best
*** Solution reached

Joining with `by = join_by(city_id)`

neotropic_communities

scree_plot(neotropic_communities)

neotropic_clusters <- kmeans(neotropic_communities[,c('NMDS1', 'NMDS2')], centers = 3, nstart = 20)
neotropic_communities$cluster = as.factor(neotropic_clusters$cluster)

ggplot(neotropic_communities, aes(x = NMDS1, y = NMDS2, colour = cluster)) + geom_point()

neotropic_cluster_hulls = convex_hull_by_cluster(neotropic_communities)

plot_nmds(neotropic_communities, neotropic_cluster_hulls, 'Neotropic', 'High')

Joining with `by = join_by(city_id)`

plot_nmds(neotropic_communities, neotropic_cluster_hulls, 'Neotropic', 'Medium')

Joining with `by = join_by(city_id)`

plot_nmds(neotropic_communities, neotropic_cluster_hulls, 'Neotropic', 'Low')

Joining with `by = join_by(city_id)`

species_in_neotropic = neotropic_communities %>% dplyr::select(city_id, cluster) %>% left_join(communities_with_traits)

Joining with `by = join_by(city_id)`

plot_phylo('Neotropic', species_in_neotropic %>% filter(cluster == 1))

resolve_neotropic = resolve %>% filter(REALM == 'Neotropic')

ggplot() +
  geom_sf(data = resolve_neotropic, aes(geometry = geometry), alpha = 0.2) +
  geom_sf(data = neotropic_communities %>% filter(cluster == 1) %>% left_join(city_points_med), aes(color = mntd_normalised, geometry = geometry)) +
  normalised_colours_scale + theme_bw()

Joining with `by = join_by(city_id, city_nm, geometry)`

plot_phylo('Neotropic', species_in_neotropic %>% filter(cluster == 2))

ggplot() +
  geom_sf(data = resolve_neotropic, aes(geometry = geometry), alpha = 0.2) +
  geom_sf(data = neotropic_communities %>% filter(cluster == 2) %>% left_join(city_points_med), aes(color = mntd_normalised, geometry = geometry)) +
  normalised_colours_scale + theme_bw()

Joining with `by = join_by(city_id, city_nm, geometry)`

plot_phylo('Neotropic', species_in_neotropic %>% filter(cluster == 3))

ggplot() +
  geom_sf(data = resolve_neotropic, aes(geometry = geometry), alpha = 0.2) +
  geom_sf(data = neotropic_communities %>% filter(cluster == 3) %>% left_join(city_points_med), aes(color = mntd_normalised, geometry = geometry)) +
  normalised_colours_scale + theme_bw()

Joining with `by = join_by(city_id, city_nm, geometry)`

Palearctic

Afrotropic

Indomalayan

indomalayan_communities = community_nmds(communities_with_traits %>% filter(core_realm == 'Indomalayan'))    

Run 0 stress 0.1193009 
Run 1 stress 0.1285877 
Run 2 stress 0.1472786 
Run 3 stress 0.1289009 
Run 4 stress 0.1240017 
Run 5 stress 0.157254 
Run 6 stress 0.1509529 
Run 7 stress 0.161852 
Run 8 stress 0.1558148 
Run 9 stress 0.1239373 
Run 10 stress 0.1165214 
... New best solution
... Procrustes: rmse 0.03581633  max resid 0.2499892 
Run 11 stress 0.1151087 
... New best solution
... Procrustes: rmse 0.02627242  max resid 0.2457618 
Run 12 stress 0.1167036 
Run 13 stress 0.1220414 
Run 14 stress 0.1229737 
Run 15 stress 0.1382712 
Run 16 stress 0.1472537 
Run 17 stress 0.1502501 
Run 18 stress 0.1666794 
Run 19 stress 0.1398437 
Run 20 stress 0.1448979 
Run 21 stress 0.1252406 
Run 22 stress 0.1443181 
Run 23 stress 0.1156951 
Run 24 stress 0.1618047 
Run 25 stress 0.1512284 
Run 26 stress 0.1245209 
Run 27 stress 0.1166941 
Run 28 stress 0.1156757 
Run 29 stress 0.1257237 
Run 30 stress 0.1169841 
*** No convergence -- monoMDS stopping criteria:
     1: no. of iterations >= maxit
    26: stress ratio > sratmax
     3: scale factor of the gradient < sfgrmin

Joining with `by = join_by(city_id)`

indomalayan_communities

scree_plot(indomalayan_communities)

indomalayan_clusters <- kmeans(indomalayan_communities[,c('NMDS1', 'NMDS2')], centers = 4, nstart = 20)
indomalayan_communities$cluster = as.factor(indomalayan_clusters$cluster)

ggplot(indomalayan_communities, aes(x = NMDS1, y = NMDS2, colour = cluster)) + geom_point()

indomalayan_cluster_hulls = convex_hull_by_cluster(indomalayan_communities)

plot_nmds(indomalayan_communities, indomalayan_cluster_hulls, 'Indomalayan', 'High')

Joining with `by = join_by(city_id)`

plot_nmds(indomalayan_communities, indomalayan_cluster_hulls, 'Indomalayan', 'Medium')

Joining with `by = join_by(city_id)`

plot_nmds(indomalayan_communities, indomalayan_cluster_hulls, 'Indomalayan', 'Low')

Joining with `by = join_by(city_id)`

species_in_indomalayan = indomalayan_communities %>% dplyr::select(city_id, cluster) %>% left_join(communities_with_traits)

Joining with `by = join_by(city_id)`

plot_phylo('Indomalayan', species_in_indomalayan %>% filter(cluster == 1))

resolve_indomalayan = resolve %>% filter(REALM == 'Indomalayan')

ggplot() +
  geom_sf(data = resolve_indomalayan, aes(geometry = geometry), alpha = 0.2) +
  geom_sf(data = indomalayan_communities %>% filter(cluster == 1) %>% left_join(city_points_med), aes(color = mntd_normalised, geometry = geometry)) +
  normalised_colours_scale + theme_bw()

Joining with `by = join_by(city_id, city_nm, geometry)`

ggplot() +
  geom_sf(data = resolve_indomalayan, aes(geometry = geometry), alpha = 0.2) +
  geom_sf(data = indomalayan_communities %>% filter(cluster == 2) %>% left_join(city_points_med), aes(color = mntd_normalised, geometry = geometry)) +
  normalised_colours_scale + theme_bw()

Joining with `by = join_by(city_id, city_nm, geometry)`

plot_phylo('Indomalayan', species_in_indomalayan %>% filter(cluster == 3))

ggplot() +
  geom_sf(data = resolve_indomalayan, aes(geometry = geometry), alpha = 0.2) +
  geom_sf(data = indomalayan_communities %>% filter(cluster == 3) %>% left_join(city_points_med), aes(color = mntd_normalised, geometry = geometry)) +
  normalised_colours_scale + theme_bw()

Joining with `by = join_by(city_id, city_nm, geometry)`

plot_phylo('Indomalayan', species_in_indomalayan %>% filter(cluster == 4))

ggplot() +
  geom_sf(data = resolve_indomalayan, aes(geometry = geometry), alpha = 0.2) +
  geom_sf(data = indomalayan_communities %>% filter(cluster == 4) %>% left_join(city_points_med), aes(color = mntd_normalised, geometry = geometry)) +
  normalised_colours_scale + theme_bw()

Joining with `by = join_by(city_id, city_nm, geometry)`

Australasia